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Abstract— This paper presents a general-purpose computer 
program which is capable of designing a large class of optimum 
(in the minim by sense) FIR linear phase digital filters. The 
program has options for designing such standard filters as low- 
pass, high-pass, bandpass, and bandstop filters, as well as 
multipassband-stopband filters, differentiators, and Hubert 
transformers. The program can also be used to design filters 
which approximate arbitrary frequency specifications which 
are provided by the user. The program is written in Fortran, 
and is carefully documented both by comments and by de- 
tailed flowcharts. The filter design algorithm is shown to be 
exceedingly efficient, e.g., it is capable of designing a filter 
with a 100-point impulse response in about 20 s. 



I. Introduction 

This paper presents a general algorithm for the de- 
sign of a large class of finite impulse response (FIR) 
linear phase digital filters. Emphasis is placed on a 
description of how the algorithm works, and several 
examples are included which illustrate specific ap- 
plications. A unified treatment of the theory behind 
this approach is available in [1] . 

The algorithm uses the Remez exchange method 
[2], [3] to design filters with minimum weighted 
Chebyshev error in approximating a desired ideal fre- 
quency response D( /). Several authors have studied 
the FIR design problem for special filter types using 
several different algorithms [4] -[13] . The advantage 
of the present approach is that it combines the speed 
of the Remez procedure with a capability for design- 
ing a large class of general filter types. While the algo- 
rithm to be described has a special section for the 
more common filter types (e.g., bandpass filters with 
multiple hands, Hilbert transform filters, and differ- 
entiators), an arbitrary frequency response can also 
be approximated. 
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II. Formulation of the Approximation Problem 

The frequency response of an FIR digital filter with 
an JV-point impulse response {h(k)} is the^-transform 
of the sequence evaluated on the unit circle, Le., 

H(fY = H(z)\ ennt = £\ h(k)e-** hf . (1) 

The frequency response of a linear phase filter can be 
written as 

H(n = G(/)e'(^-(^)2*f) (2) 

where G(f) is a real valued function and L = 0 or 1. 
It is possible to show that there are exactly four cases 
of linear phase FIR filters to consider [1] . These 
four cases differ in the length of the impulse response 
(even or odd) and the symmetry of the impulse re- 
sponse [positive (L = 0) or negative (L = 1)] . By 
positive symmetry we mean ft(fe) = h(N - 1 - ft), and 
by negative symmetry h(k) = - h{N - 1 - ft). 

In all cases, the real function G(f) will be used to 
approximate the desired ideal magnitude specifica- 
tions since the linear phase term in (2) has no effect 
on the magnitude response of the filter. The form of 
G{f) depends on which of the four cases is being 
used. Using the appropriate symmetry relations, 
G(f) can be expressed as follows. 

Case 1: Positive symmetry, odd length: 

G{f) = Z a(k)cos(2nkf) (3) 
ft = o 

where n = (N- l)/2, o(0) = h(n), and a(k) = 2h{n - ft) 
for ft = 1, 2, • • , n. 
Case 2: Positive symmetry, even length: 

G(fl= t b(k)cos[2n(k~±)f] (4) 

where n= N/2 and b(k) = 2h(n - fe)forft = l, 
Case 3: Negative symmetry, odd length: 

Gif) - L <*k) sin (2*fc/) (5) 
ft = i 

where n = (N - l)/2 and c(ft) = 2ft(n - ft) for ft - 1, 
2, * • , n and h{n) = 0. 
Case 4: Negative symmetry, even length: 

Gif) - £ [2ir(fe - -|)fl (6) 

where n = N/2 and d(ft) = 2h(n - ft) for ft = 1, • • • , n. 

Earlier efforts at designing FIR filters concentrated 
on Case 1 designs, but it is now possible to combine 

1 For convenience, throughout this paper the notation H(f) 
rather than Hie* 2 *?) is used to denote the frequency response 
of the digital filter. 
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all: four cases into one algorithm. This is accom- 
plished by noting that G(f) can be rewritten asG(f) = 
Q(f)P(f) where P(f) is a linear combination of cosine 
functions. Thus, results that have been worked out 
for Case 1 can be applied to the other three cases as 
weli. For these purposes, it is convenient to express 
the summations in (4)-(6) as a sum of cosines di- 
rectly. Simple manipulations of (4)-(6) yield the 
expressions. 
Case 2: 



£ 6(fe)cos[2tr(ft--|)f] 

-coB(»/)j| 6(fe)cos(2irfef). (7) 

Cose 3: 

T c{k) sin (2**/) - sin (2rf) £ c(k) cos (2trfcf ). 



Case 2: 4 



(8) 

Case 4: 
± d(fe)sin[2tr(fc--|)n 

k = l 

= sin (nf) £ d (k) cos (2nkf ) (9) 

k = 0 

where 

r b(i) = 6(0)+-| &(D 
6(*)--|[S(*- + 

fe = 2, 3, • • • , n - 1 
J>{n)=\Hn- 1) (10) 

f c(l) = £(0)--|c(2) 

C(fe)=i[0(fe- 1)-C(* + 1)], 

= 2, 3, • • • , n - 2 

c(n- l)=lc(n- 2) 

c(n) = \Z(h- 1) (11) 

fd(l) = d(0)--|d(l) 

d(ft) = |[S(fe- l) -5(*)j, 

fe»2, 3, • • ,rt- 1 

k d(rt)-id(/i- 1). (12) 

The motivation for rewriting the four cases in a 
common form is that a single central computation 
routine (based on the Remez exchange method) can 
be used to calculate the best approximation in each 
of the four cases. This is accomplished by modifying 
both the desired magnitude function and the weight- 



Case 3: < 



Case 4: < 



ing, function to formulate a new equivalent approxi- 
mation problem. 

The original approximation problem can be stated 
as follows: given a desired magnitude response D(f) 
and a positive weight function W(f), both continuous 
on a compact subset F C [0, \ ] (note that the sam- 
pling rate is 1.0) and one of the four cases of linear 
phase filters [i.e., the forms of G(f)] , then one wishes 
to minimize the maximum absolute weighted error, 
defined as 

- max W(f) \D(f) - G(f)\ (13) 

over the set of coefficients of G( f )• 

The error function E(f) can be rewritten in the 
form 



Eif) - W(f) [D(f) - G(f)} = W(f) Q{f) 



-p(f)\ 

(14) 



if one is careful to omit those endpoint(8) where 
Q(f) = 0. Letting £(/) = D(f)/Q(f) and W(f) = 
W(f) Q(f), then an equivalent approximation problem 
would be to minimize tile quantity 

\\E(f)\\ = max %(f ) \3(f) - P(f)\ (15) 

by choice of the coefficients of P(f). The set F is re- 
placed by F' » F - {endpoints where Q(f) = 0}. 

The net effect of this reformulation of the problem 
is a unification of the four cases of linear phase FIR 
filters from the point of view of the approximation 
problem. Furthermore, (15) provides a simplified 
viewpoint from which it is easy to see the necessary 
and sufficient conditions which are satisfied by the 
best approximation. Finally, (15) shows how to 
calculate this best approximation using an algorithm 
which can do only cosine approximations. The set 
of necessary and sufficient conditions for this best 
approximation is given in the following alternation 
theorem [2] . _ 

Alternation theorem: UP(f) is a linear combi- 
nation of r cosine functions i.e., 

P(f) - 2 a(k) cos 2irkf, 

then a necessary and sufficient condition that 
P(f) be the unique best weighted Chebyshev ap- 
proximation to a continuous function D(f) on 
F' is that the weighted error function E(f) = 
$(f) 0(f) - P(f)] exhibit at least r + 1 extre- 
mal frequencies in F' . 

These extremal frequencies are a set of points 
{F t }, i - 1, 2, • • • , r + 1 such that F x < F 2 < • • • < 
F r < F r + X , with E(Fi) * - E(F Ul ), i « 1, 2, • - ■ , r and 

An algorithm can now be designed to make the 
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INPUT FILTER 
SPECS 



-FILTER 



[BANDPASS 



TYPE ■ 



1 differentiator! 



WATE 
WEIGHT FUNCTION 



HILBERT 
TRANSFORM 



DENSE GRID FOR 
THE SET F 



EFF 

DESIRED MAGNITUDE 



FORMULATE EQUIVALENT 
APPROXIMATION PROBLEM 



INITIAL GUESS FOR r+1 
EXTREMAL FREQUENCIES 



REMEZ 
SOLVE THE 
APPROXIMATION PROBLEM 



CALCULATE IMPULSE 
RESPONSE 



PRINT OUT THE 
OPTIMAL ERROR a 
IMPULSE RESPONSE 



Fig. 1. Overall flowchart of filter design algorithm. 



error function of the filter satisfy the set of necessary 
and sufficient conditions for optimality as stated in 
the alternation theorem. The next section describes 
such an algorithm along with details as to its imple- 
mentation. 

III. Description of the Design Algorithm 

As seen in Fig. 1, the design algorithm consists of 
an input section, formulation of the appropriate 
equivalent approximation problem, solution of the 
approximation problem using the Remez exchange 
method, and calculation of the filter impulse re- 
sponse. The flowcharts of Figs. 2-5 give details of 
the exact structure of the computer program. 

The input which describes the filter specifications 
consists of the following. 

1) The filter length, 3 < nfilt <nfmax (the upper 
limit set by the programmer). 

2) The type of filter (jtype): 

a) Multiple passband/stopband ( jtype =1) 

b) Differentiator (jtype =2) 

c) Hflbert transformer (jtype =3). 

3) The frequency bands, specified by upper and 
lower cutoff frequencies (edge array) up to a maxi- 
mum of 10 bands. 

4) The desired frequency response (fx array) in 
each band. 

5) A positive weight function (wtx array) in each 
band. 



6) The grid density (lgrid), assumed to be 16 unless 
specified otherwise. 

7) Impulse response punch option (jpxjnch). 

Part 3) specifies the set jF to be of the form F = UB, 
where each frequency band B t is a closed subinterval of 
[0, The inputs 4) and 5) are interpreted differ- 
ently by the program tor a differentiator than for the 
other two types of filters (see the eff and wate sub- 
routines in Figs. 3 and 4). The weight specification in 
the case of a differentiator results in a relative error 
tolerance as is used in all other cases. 

The set F must be replaced by a finite set of points 
for implementation on a computer. A dense grid of 
points is used with the spacing between points being 
0.5/(lgrid X r) where r is the number of cosine basis 
functions. Both D(f ) and W(f) are evaluated on this 
grid by the subroutines eff and wate, respectively. 
Then the auxiliary approximation problem is set up 
by forming D(f) and $(f) as above, and an initial 
guess of the extremal frequencies is made by taking 
r+1 equally spaced frequency values. The subrou- 
tine remez (Fig. 5) is called to perform the calcula- 
tion of the best approximation for the equivalent 
problem. The mechanics of the Remez algorithm will 
not be discussed here since they are treated elsewhere 
tor the particular case of low-pass filters [9] . (The 
flowchart of Fig. 5 gives details about the mechanics 
of the Remez algorithm as implemented in this design 
program.) 
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input fitter specs: 

1. Filter Length (NFILT) 

2. Type of Filter (JTYPE) 

3. Number of Frequency Bands (NBANDS) 

4. cord Punch Option (jpunch) 
sl orid Density (low) 





Input filter Specs : 

1. Bondedges — Array E0GE{ ) 

2. oeelred Function or Desired slope 
(Differentiator) In each bond-Array FX{ ) 

3. weighting co nst a nt for eoch band— Array 
WTX( ) 



i 



WATE to 
calc ul ate WT 



Calculate the Number of Approximating 
Functions (MFCpJS). This depends on the 
filter type and filter length 



Set up o dense grid array ORID with spacing 
between grid points eqool to 
aS/dJBRlDXNFCNS) 



FX array 



Set up two arrays OES< J and WT( ) 
containing the desired function Mhie and 
the weight function at each grid point 



Type i (etuttibandl 



Odd Length 



1 



SUBROUTINE 
CFF to 
caScutote OCS 



* Types 2ondS(i 



Even Length 



test i ) -de* j Vccef » «rbx i >J 
wT(j)=wT(J) cosf^eRiw j >] 



.Even Length 



OES(jl»DESU)/pn[ vscm))] 

wTij)=wT<j) sin[ Texan])] 



. Odd Length 



0E9jHXS(j>/sin[ ZTQROXli 
WT(j>VTIj) sin[ zw&mxo] 



Moke an initial guess for the location of 
the extremal frequencies. The location is 
recorded by keeping track of the index of the 
frequency in the 6Rto array. 
The indices ore stored in IEXT array. 
Initial guess is nfcns + 1 equally spaced 
index values 



I 

Fig. 2. Detailed flowchart for fitter design algorithm. 



The appropriate equations (3)~(12) are used to re- 
cover the impulse response from the coefficients of 
the best cosine approximation obtained in the bemez 
subroutine. Hie outputs of the program are the im- 
pulse response, the optimal error (min \\E(f)\\), and 
ther + 1 extremal frequencies whereof) = ±||£(f )||. 

It is possible that one might want to design a filter 
to approximate a magnitude specification which is 
not included m the scheme given above, or change the 



weight function to get a desired tolerance scheme. A 
flowchart of such a program is given in Fig. 6. In 
such cases, the user must code the subroutines eff and 
wate to calculate D{f) and W(f). The input is the 
same as before, except that there are only two types 
of filters, depending on whether the impulse sym- 
metry is positive or negative. 

A detailed program listing of the generalized design 
program is given in the Appendix. Representative 
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1 



Coll REMEZ exchonge algorithm 



Type 1 



Odd Length 



H{ j)= 1/2 Ct(NFCNS + 1-j) 
H(NFCNS) = a(1) 
Where a array contains 
coefficients of best 
cosine approximation 



Type 2 or Type 3 



Even Length 



H(j) *1/4 [a<NFCNS + 1-j> + 

a( NFCNS + 2 — j )] 
HID =1/4 CMNFCNS) 
H(NFCNS> = l/2a<1) + 
1/4Q(21 



Even Length. 



H(j) =1/4[a(NFCNS + t-j)- 

a( NFCNS +2-j)] 
H(1) =1/4 a(NFCNS) 
H(NFCNS)M/2a<!)- 
!>4 a(2) 



Odd Length 



H< j)=1/4 [a(NFCNS+1-j)- 

a( NFCNS +3-j)] 
H(1)=1/4 O(NFCNS) 
H(2)«l/4a(NFCNS-i) 
H ( NFCNS ) * 1/2 Ct ( t) - 

1/4Q(3) 
H (NFCNS + t) =0 



Program Output 

1. Title 

2. Type of Filter 

3. Filter Length 

4. impulse Response 

5. Bandedges 

6. Desired value/ slope In each band 

7. weighting in each band 

8. Deviation In each band 

9. Deviation in ds for Type i 
to. Extremal Frequencies 



Punch impulse response if 
desired, JPUNCHH 



Fig. 2. (Continued.) 



SUBROUTINE EFF 




Evaluates desired function at o grid point 

Fig. 3. Flowchart fox subroutine eff. 
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SUBROUTINE WATE 



t OR 3 




WATE(F) * WTX(J) 
for F c jth bond 



WATE(F) 



WTX(j) if FX<j>«0 
WTXiJVF if FX(j)*0 



T 



Evaluotes weight function ot o grid point 

Fig. 4. Flowchart for subroutine WATE. 



REMEZ 



(tS) Main 



Set maximum number of iter o lions - 

ITRMAX-29 
Initial tee iteration count 

NITER «0 
Let n * • of approxifftottng functions 



Iteration Loop -Statement too 



NITER* NITER + 1 






60 to 400 




calculation of 




coefficients of 




best approximation 



Calculate the abscissae for the 
Logronge tntcrpotatton 
X(j)>O0s[2srF(j)] 

Where Ft}) ore the extremal 
fregjusncies 




r 


calculate Lagrange 
coefficients usfng i 


> Interpolation 
Subroutine D 

~X(j)] (i#j) 



Calculate the current value of the deviation 
(OEV) and fl prlnt DEV. 

£ AO(J) D£S[XEXT0)] 

oev- -i= « 

I f-t)J"»«X j l/wr[lEXTt i i] 



Hi 



[ Record son CDEV) 


»setO£V»|D£V| 1 




> 


Calculate ordlnctes tor Lagrange interpolation 

y(i)-OEs[IEXT<J)]♦HI , DEV/WT[lEXTCi)] 






OUCH 
error 
message 






: •> 



60 to calculotion 
y of coefficients 
of best 

approximation 



Fig. 5. Detailed flowchart for subroutine REMEZ. 
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For the tst extremal frequency set upper 
kup otk> lower KLOW limits on the 
Grid points to be searched 

KUP -XEXTC2) 

KLOW-0 

Furthermore, let the sign of The error ol 
each extremal frequency be denoted 

by<rti) 






using 


Subroutine GEE at GfttD <k+11 


where 


k-iexti i } (t.a. me erto poi 


nt arfjocent 


to the ]th extremal ffeouenc] 


rk Ttoeo 


coteulate the weighted error 


B{K+1)] 


ERB-WTtK-HlfeEEtK+IJ-DE 



ree Oo to endpoint search [Statement 300) . 



AO X V . arrays 



1 i i 



SUBRQVTtttC €EE 



2 AO tup YW 

_ n&t x-xttn) 




JL, x-xtml 

Where 

. x»cos[2*-grio(ki] 



w fl ighte d error at index it— I 



There is a toca 


1 moijnwni of 


the error cone 


where the 


signed error is 


■water thort the 


present deviation. Continue 


the search at K+3, 


...KUP ontil this local max 


is found. 





Them fs o local ma* of 
error carve, search K-2. 

K-9 KLOW until tMs toed 

max is found 



1. Change this extremal frequency 

2. update KLOW and kup for 



(tnewlExn 
KUP.min{lEXT(Jf2), NGRto} 
5. note that o 



Oo to beginning 

of the bop 
[Statement £00) 




conffnee eeorcNna k-2,k-s, 
— KLOW for a local max 
error is 

tPSV 




Search *+2.K+S,— kup far 
a taca) max where signed 
error tOEV 



^oon't find 
«o to £00 



Fig. 5. (Continued.) 
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Statement 300-Endpotnts March 



Oont find such o 
loco) max 



1 



search those indices less than mm {iEXT(i), new IEXTO )} 
for error with sign ~o*(D and greater than the error 
at new lEXT(n) 



Find such a 
local max 



| Store the index and error value | 



— „ ., W IM I— — III. I I —.■„... ..,„ 

Search those indices greater than mo*{lEXT<nl,n«w lEXTtn)} 
for local max of the error curve with sign -<r(n) 
and error t error at new r EXT CO 
and error found above 



Don't find such a 
max 



I EXT for fw 


act 


iteration 


IEXT (D-d 


*h 


ex found above 


IEXT (2) « t 

i 


ten 


t iexto) 


IEXT(n) =r 


tew 


IEXT (n-D 


1 


r, . 




Find such a 
local max 



| IEXT Ql* new IEXT l)\ I 



IEXT (I) -new IEXT (2} 

IEXT (OH) - new IEXT (n) 
IEXT(nl » index just found 



t 



L 




Yes 



-»>6o to 100 



statement 400 



Evaluate SEE at KFCNS equally 
points in the interval [o.t/2] 
f n l n-t \ 



Inverse DFT to obtain coefficients of the 
best cosine approximation 



Fig. 5. (Continued.) 



input Fitter specs - 

NFILT, JTYPE.NBANDS, JPUNCH. LGR1D, 

EDGE ( ), FX( ), WTXJ ) 



FX array 



1 



SUBROUTINE 
EFF 

Supplied by 
user 



1 

calculate number of a 
and set up the dense < 


r„ — 1 

pproxi mating functions 
irtd 




L . 


set up the arrays 0€ 


9 and WT 



Positive Symmetry 
( Type I ) 
— * 



Odd Length 



Negative Symmetry 

(Type 2) 
P 



Even Length 



Odd Length 



continued as in the main flowchart 

Fig. 6. Flowchart for arbitrary magnitude filter 

ritfam. 



wTX orroy 



1 



SUBROUTINE 
WATE 

Supplied by 
User 



Even Length 



design algo- 
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LINEAR PrtASt DIGITAL FIcTtft UcSXGN 
RztttZ cAurtAh lit ALGORITHM 
aANOPASo flULR 



***** IM»»UL3t K£SPJWS£ 



H( 
HI 
H( 
H( 
H( 
Hi 
H< 
HI 
HI 



1) 
2) 
3) 
*) 
$1 
6) 
7) 
3> 

9) 



0.337*b9«7c> 
0.1*»93*>299t- 
0.1U569360t- 
0.25<il50&7c- 
-0.15929J92c- 



H( 10) 
H( 111 
HI 12) 



02 
01 
01 
02 
01 



HI 
HI 
HI 
HI 
HI 



-0.3<»0a53<t3c-01 = tyl 



-Q.38112i77t-01 
-0.1*625i69c-01 
a.*00tt9ti<*lt-01 

a«u*i»oriiiL ao 

0.1«id$Q752t 30 
0.233S**o06c U0 



HI 
HI 
HI 
HI 
tit 
HI 



£-»> 
23) 
l£> 
21) 
20) 
19) 
ii) 
17) 
Ibi 
19) 
j.<») 
li) 



LOHcK dANU cOGc 
UPPcR dANU ctGc 
DfcSIfcLO VALUE 
HcIGHTIrfG 
UEVIATION 



BAND x 

0. 

u. OdocyouG 

±. QU30UGJ0 

i* oooooooo 

0. 012*33oh 



QANO 2 
O.lbO iOOGfl 

d.5ti jauaou 
j • 

x.QOuJOOOJ 



dANU 



UtVlATIJN Itf DB -3d. 1 Jo 03^+ 13 -3d. 10d03*»13 



€*TtitHAL F*cGU£NCitS 

0. &.036<»583 0.0o770d3 0,0600000 3.1600000 

G. 1730206 0.2Qod750 Q.2*?937? U.^37o0<*2 J.3<>idS5Q 

0.3737500 0.4256251 G*J>75104»3 



TIME 3 0.7&9MI63 SECONDS 

Fig, 7. Output listing for an N= 24 low-pass filter. 



UOWPASS FILTER 
N*24 




0.2 03 
FREQUENCY 



Fig. 8. Magnitude responses, on linear and log scales, for an 
N = 24 low-pass filter. 



input card sequences are given for the design of a 
bandpass filter and a differentiator. To approximate 
an arbitrary magnitude response and/or an arbitrary 
weighting function, all the user has to do is change 
the subroutines epf and wate and use the program in 
the Appendix. In the next section, representative 
filters designed using these algorithms are presented. 

IV. Design Examples 

Figs. 7-22 show specific examples of use of the de- 
sign program for several typical filters of interest. For 
each of these filters, one figure shows the computer 
output listing (including the run time on a Honeywell 
6000 computer), and the other figure shows a plot of 
the filter frequency response on either a linear or a 
log magnitude scale (or sometimes both). Figs. 7 and 
8 are for an N = 24 low-pass filter. For this example, 
the run time was 0.77 s. Figs. 9 and 10 are for an 
N = 32 bandpass filter. This example is the first ex- 
ample listed in the prologue to the program in the 
Appendix. The run time for this example was 0.82 s. 
Figs. 11 and 12 are for an N - 50 bandpass filter in 
which unequal weighting was used in the two stop- 
bands. Thus the peak error in the upper stopband is 
ten times smaller than the peak error in the lower 
stopband. A total of 2.96 s was required to design 
this filter. Figs. 13 and 14 are for an N = 31 bandstop 
filter with equal weighting in both passbands. For the 
design of this filter 1.61 s were required. 

To illustrate the multiband capability of the pro- 
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FINITE INPULSc RESPONSE CFI*) 
LINEAK PHASE DIGITAL FIlT£R DESIGN 
RtNtZ EXCHANSc ALGCWITHN 
BANOPASS FILTER 
FILTER LENGTH = 32 

IHPULSE RESPONSE ♦♦♦♦♦ 



LOWER BAND EOGE 
UPPER BAND £064 
OESIREO VALUE 
WEIGHTING 
OEVIATIOM 



H( 


11 


s 


-0*57:*34Ai:i£-02 


8 


H( 


32) 


H( 


2) 




0.49027198E-03 


S 


HI 


31) 


H( 


31 


* 


4*757o3S«»5E-02 


8 


HI 


30) 


H( 


4) 




*0.65l4ll92£^02 


8 


H< 


29) 


H( 


5) 


a 


0.13960525tH)l 


8 


H( 


28) 


R( 


6) 


«. 


(U22951469E-02 


8 


H( 


27) 


H( 


7» 


3 * 


-0*1 99940 67s.- 01 


m 


H« 


26> 


HI 


8) 




Q.7M6956QE-02 


8 


HI 


2S) 


H( 


9) 


8 


-Q.39657363E-01 


8 


Hi 


24) 


Hi 


10) 


s 


0.1126Q114£-Qi 


S 


Hi 


23) 


H( 


11) 




0 • 662336 *3E- 01 


a 


Hi 


22) 


H( 


12> 




-b;iO^97223E*01 


s 


Hi 


21) 


H< 


13) 


a 


0*85J.36133£-0i 


s 


Hi 


20) 


H( 


14) 


s 


-0.12024993E 00 


* 


H( 


19) 


H( 


15) 


s 


-0*29676577E 00 


8 


Ht 


IB) 


H( 


16) 


a 


0.30410917E 00 


8 


H< 


17J 




BAND 


1 BANO 2 







0*10000000 

0* 

10*00000000 
0* 001513*2 



DEVIATION IN 06 -56*40254641 



0*20000000 
0*35000000 

1*600000*0 
i,oooaoopa 

0« 01513110 



BANO 3 
0*42500000 
0*50000000 
0* 

io: 00000000 
0*00151312 



BANO 



-36*40254641 -56.00254641 



EXTREMAL FREQUENCIES 

0* 0*0273437 0*0527344 

0*1000000 0.2000000 0*2195312 

0*3132612 0*3386719 0*3500000 
0*4503906 0. 4796675 



0*0761719 0*093750p 
0*2527344 0*2839644 
0*4250000 0*4320125 



TINE* 0*6245625 SECONpS 

Fig. 9. Output listing for an itf = 32 bandpass filter. 



BANDPASS FILTER 
W»32 




0.2 OS 
FREQUENCY 



Fig. 10. Log magnfode response for an N= 32 bandpass filter. 
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FINIT£ iHPUtit *tSKi*Sfc 

uiucAR PHASt OIGITAu FXlT£R OcSiS. 1 * 

<£MtZ tXCriAUoc ALiiUKITHH 

34 NOP A £S FitTs.* 
F iLTclR LcNSrn = 55 
• XtlPUt-St. RcSPQiJSt 



Hi 11= J»lJ~Ob£3>£C-0£ 9 fit 

M< 2* = a,637776i.5t-0<; * H( 

H( 3) = Q,3$?55olQc-02 - HC 

HI 4) = -0. j0b77i55t-02 = HI 

H( 51 = -0. 93*063796-02 = H( 

H( 6) = 3.2Jl55o29t-U2 * HI 

H< 7) « D.3*o37965=>02 a H< 

hi a* = o,iii7zasoc-oi * hi 

HI 9) = 0.1tb467S9c-01 * HI 

HI 10) « -0.39o307<Ji»t-02 * HI 

H( Hi = -0,923a4245c-02 * HI 

HI 12) = -0.20*»06392c-01 = HI 

HI 131 = -G.19<f6Q463c-61 » HI 

HI 14) « 0.312430l3t-0i * HI 

HI 15) = O.b3Q455b7£-02 » Hi 

HI 16) = -0.20462303E-01 * HI 

HI 17) » O.o5740513c-02 - HI 

HI 16) = -0.11202i27t-Q2 « H( 

HI 19) = 0.41956965£-01 « HI 

HI 20) = 0.35764266£-01 = HI 

HI 21) = 0.34744602E-Q1 * HI 

HI 22) = 0.71496359t-01 = HI 

HI 23) = -Q.17138631£ U0 * Hi 

HI 24> = -0.ia235d44£ 00 = H( 

HI 251 = 0.74Q5S02*£-01 = HI 

HI 26) = -0. 10317421c d0 = HI 

HI 271 = a.2S716721fc-0l * Hi 

HI 281 = 0.376135*fr7fc OO = HI 



>5) 
5*) 
53) 
52) 
51) 
50) 
49) 
46) 
47) 
46) 
45) 
44) 
43) 
4>) 
41) 
40) 
39) 
3d) 
371 
36) 
35> 
34) 
331 
32) 
31) 
30) 
29) 
26) 



. . 6AM D 1 

LOHcR BAND tDGE. 0* 

UPPcR BAND £OGc 0.05000000 

QtSIRtD VAlUc 0. 

WEIGHTING 10.00000000 

OlVIATIOn 0.003v4*do 

OtVlATIiiN IN 06 -*9.<:Sb57a3<> 



SAND 2 

u.ipoooooo 

U.lsCOOOGO 
a.OOOOOuOu 
1. OUoJOQ&Q 
0, 03«****»at>9 



3 AMD 3 

d.iaoudodo 

0*^5000000 
0. 

3.00000000 
il.Qii**8*do 



-iJ.25o57034 -39.79499**9 



dA>40 if 
u.30000000 
0*36000600 
1*00600000 
1*04000000 
O.Q34*h6>9 
-c9 .t5o570ji» 



LOWcK BAND tt)G£ 

UPPcR BAND cOGt 

OcSIREO VACUt 

WEIGHTING 

0&VIATION 

Ot VI AT ION IN OB 



BAND 5 
0.410000JO 
0*50000000 
0. 

20. 00000000 
0.001722*3 
-55.27717018 



dANO 



£XTRlHAl FRfcQUtNCIcS 



0*1000000 

o.iaooooo 

0.2436160 
0.350223i 
0.4457143 



0.0167*11 
0.10692 66 
0*16556 04 
0.2500000 
0*3600000 
0.4635714 



0.0323o61 
O.l2o7o57 
0.1976571 
0.3300000 
O.ifiOdOOO 
0.4614265 



6.0446429 
0 .l4it<»^07 
0*£l5462l 
0.3122766 
0.41&56O4 
0.500U00O 



0*0500000 
0.1500000 
6*2302232 
0.3323661 
0«*2d9732 



TIME= 3.8161(219 SECONDS 

Fig. 15. Output listing for an N= 55 multiband filter. 



MULTfBAND FILTER 




02. 03 
FREQUENCY 



Fig. 16. Log inagnitu.de response for an N = 55 multiband 

filter. 
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C 

C THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128 « BUT 
C THIS UPPER LIMIT CAN BE CHAMGEO BY REOINLNSIONIN? THE 
C ARRAYS lEXTt AOt ALPHA* X« Y* H TO BE NFMAX/2 ♦ 2. 
C THE. ARRAYS DES, 6JU0, AND NT MUST OXMEMSIONCO 
C 16<NFNAX/2 +2). 
C 

NFMAXsUe 
100 CONTINUE 
JTYPESO 

C 

£ PROGRAM INPUT SECTION 
C 

RE AO * * NFILT , JTYPE t NBANOS ♦ JPUNCM t LGRIO 
IFCNFILT.GT.NFNAX.01UNFILT.LT.3) CALL ERROR 
IF C N8AN0S *LE.0> NBANOSsl 

C 

C GRIO. DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED 

C OTHERWISE 

C 

IF(LGRIO.LE.O) L6RI0=16 
JBs2«MBANOS 

REAO *t(EDGE«J)*U*ltJB) 
READ *, I FXiJ)tU»i •NBANOS) 
REAO *?(WTX<«J) *«Jsl« NBANOS) 
IF(JTYPE.EQ«0) CALL ERROR 

NIEGs; 

XFtJTYPE.EQ.l) NE6eO 
NOOOsNFILT/2 
NOOOsNFI LT-2*NOOO 
NFCNS=NFlLT/2 

IFtNOO0.EQ.l«ANOiNE6.EQ*0) NFCWS=NFCNS*1 

C 

C SET UP THE DENSE GRID* THE NUMBER OF POINTS IN THE GRIO 

C IS (FILTER LENGTH ♦ l)*6RIO DENSITY/2 

C 

GRIO|l)=EOGEU) 
OELFsLGRIO*NFCNS 
OELFsO,5/OELF 
X.FtNEG.EQ.01 GO TO 133 
IFfEOGE(l).LT.OELF) 6RI0ll)=0ELF 
13S CONTINUE 
J=l 
L=l 

LBANOsl 
140 FUPbEOGEIL+1) 
145 TEHPsGRiDCJ) 

C CALCULATE THE OLSIREO MAGNITUDE RESPONSE AND THE *E*6HT 

C FUNCTION ON THE GRID 

C 

□ES ( J ) eEFF I TLRP » FX « NTX * LOANO f JTYPE > 
WT< J) ATE « TEMP »FX« WTX « LBANO* JT YPE > - 
u=J*i 

GRIOXJJ=TCMPiDtLF 
IFtGRIUtUl.GT.FUP) GO TU 150 
GO TO 145 
1SU GRIOCJ-1 )=FUP . 

DES( J-l)=E.FF(FUP.Fi.*TX.LUANO< JTTPt) 

tfY t J-l 1 ew ATE ( FUP t F X t V TX t LB AND t JT YPE ) 

LBANO=LBAND+l 
LsL*2 

IF tLBANO.GT. NBANOS ) GO TO 160 
GRlD<J)*fcOGE<L) 
GO TO 140 
160 NGRIOsJ-1 

IF(NEG.NE.NODO) GO TO 166 

IFCGRIOCNGRI0)»GT.<0.5-DELFI) NGRlUsNGRlD-1 
165 CONTINUE 

C SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT 

C TO THE ORIGINAL PROBtEM 

C 

IF(NE6I I70tl70«160 
170 IF(NODO.EO.l) GO TO 200 

DO 17S J^ltNGRIO 

CHAN6Es0C0S<PI*6RI0(U> 1 

OES (U)=0ES(«J) /CHANGE 
175 tfT<J)s*T(«J)*CHAN6£ 

GO TO 200 
180 IFtN00O.EO.il GO TO 190 

DO 105 UsltNGRIO 

CHAN6EsQSIN(Pi»GRIDCU> ) 

OES < U ) sOES f U ) /CHANGE 
165 tfT(dl=MT<J)*CHAN6£ 

GO TO 200 
190 UO 195 J=l.NGRIO 

CHANGEsOSlN<PI2*GRIDIU» I 

DES I U > sOES I J I /CHANGE 
195 VTf J)£|fT<J)*CHANGE 

C 

C INITIAL GUESS FOR THE EXTREMAL FREQUENCIES— EQUALLY 

C SPACED ALONG THE GRID 

C 

200 TENpsFLOAl (NGR ID*I ) /FLOAT (NFCNS J 

00 210 NFCNS 
210 IEXT<J)s<J-ll*TEMP+l 

IEXT C NFCNS*! 1 sNGRIO 

NNlsNFCNS-1 

N2=NFChS*l 



C 

C CALL THE KENEZ EXCHANGE ALGORIlHM TO 00 THE APPROXIMATION 

C PROBLEM 

C 

CALL REME2 1 EU6E * NBANOS) 

C 

C CALCULATE THE IMPULSE RESPONSE. 
C 

IFINEG) 300«<*00.320 
300 IFINOOU.EU.Ol GO TO ilO 

UO 305 J=1»NM1 
605 H4JlsO«5*ALPHA(N2-u> 

HlNFCNSIsALPnA(l) 

GO 10 350 
310 H( 1 ) =0.25* ALPHA tdlFCNli I 

UO 31 5 Js2«WU 
ili H<J)=0.25*iALPHA<N*-U)+ALPHA<NFCNS*2-0ll 

»l( Hf CUH l su . 5* ALPh A C 1 ) 1 0 . 25* ALPHA ( i ) 

GO TO 350 
320 IFINuUu.Efe.O) GO 10 33u 

Hll)=U.25*ALPHA<NFCNb> 

H < 2 1 =0 »25»ALPnA « NM1 > 

00 325 J=3*Nhl 

329 H( J)s0.25*(ALPHA|Mic-J»-ALPHA(NFCNS43-JII 
H(NFCtt»)=0*S*ALPHAU>*0«2S*ALPHA(31 
H(NZ»=0.0 

GO TO 350 

330 HCl)sQ»2S*ALPHAlNFCNS> 
DO 335 J=2 tNMi 

335 H I U I so • 25* ( ALPHA t N2- J ) -ALPHA ( NFCNS*2-U) ) 
H I NFCNS ) =0 • 5* ALPHA < I J -0 ♦ 25* ALPHA < 2 > 

C 

C PROGRAM OUTPUT SECTION. 

C 

350 PRINT 360 

360 FORMAT tlHl# 70 <1H*)//25X» 'FINITE IMPULSE RESPONSE (FIRM/ 
125Xf * LINEAR PHASE OIG1TAL FILTER DESIGN*/ 
225Xi«REME2 EXCHANGE ALGORITHM*/ ) 
IFIJTYPE.EQ.il PRINT 365 
365 FORMAT 1 25X t * BANDPASS FILTER*/) 

IF<UTYPE.E0.2) PRINT 370 
370 FORMAT I 25Xt« DIFFERENTIATOR*/) 
IF(JTYPE.EQ.S) PRINT 375 

375 FORMAT ( 25 X t * HILBERT TRANSFORMER • / 1 
PRINT 376 » NFILT 

376 FORMAT ( 15X t * F ILYLR LENGTH » *,I3/) 
PRINT 360 

360 FORMAT ( 1SX t • ***** IMPULSE RESPONSE **•**•) 
00 361 Usl t NFCNS- 

K=NFILT+1«-J 

IFtNEG»E0,O> PRINT 362*UfH|J)«K 
IFCNEG.EO.l) PRINT 363tU.HlJ)tK 

361 CONTINUE 

362 FORMAT (20X? "HI * t IS« * ) -» *fE15.6,* * HI,M»t')'l 

363 FORMAT <20X**H< * *I3< *) = **E19*6t* * -Ht*»I4**)*) 
IFINEG.EO.i.ANO.NODO.EO.i) PRINT 364»Wt 

364 FORMAT(20Xt«H(*rI3*M = O.OM 
00 450 Kelt NBANOS* 4 

KUPsK>3 

IF (KUP.6T. NBANOS) KUP=N BANDS 
PRINT 36S.(J«J=K.KUP> 

365 F0RNAT</2<tXt4i»BAMD'tI3«6X)> 
PRINT 390*<ED6E<2*U-i)fUsK,KUP> 

390 FORMAT <2Xt 'LOWER BAND EOGC* «5F15.9> 

PRINT 395«(EDGE(2*U)*J«KtKUP) 
395 FORMAT ( 2X t * UPPER BAND EDGE* *5F15»9t 

IFIUTYPE.NE.2) PRINT 400. <FX(U) , JsKfKUP) 
400. FORMAT C 2X * * OESIREO VALUE • ? 2X 1 5Fl5 »9) 

IFCJTYPE.E0.2) PRINT 405, IFXI J) t J=K.K UP I 
403 FORMAT 1 2X * * DESIRED SLOPE* «2Xt5Fl5t9) 

PRINT 410*(WTX(U>tU«K*KUPI 
410 F0RMAT|2X**VEtGHTlNGt»6X«5F15.9) 

DO 420 J=K.KUP 
420 0EVIAT(JJcOEV/WTXU) 

PRINT 429t(0EVlAT(J>tU3KtKUP) 

429 FORMAT ( 2X t * OEV I AT ION • # 6X * 5F15. 9) 
IFCJTYPE.NE.l) GO TO 450 

UO 450 U=K.MjH 

430 DEVlAT(J)s20.0*ALO61O(Ut.VlA1{JI ) 
PRINT <»35*<&EVlAT(U)tj5K«KuP) 

435 F0RNAT(2Xt'uLVlATl0N IN De««SF15.9) 
450 CONTINUE 

PRINT 455t CGHlO(lLXTlO))t«lsl,NZ) 
455 FORMAT i/2X,»LXTKEAAL FRLUUENCIESVC2X*5F12.7) I 

PRINT 460 
460 FORMAT(/lXf70tiH*)/lHl) 

IFCJPUhCH.NE.O) PUNCH *• (H|U) • J=l » NFCNS) 

IFINFILT.NE.O) GO TO 100 

RETURN 

ENU 

FUNCTION EFFITEMP.FXtNTX.LBANOtJTYPE) 

C FUNCTION TO CALCULATE THE OESIREO MAGNITUDE RESPONSE 

C AS A FUNCTION OF FREQUENCY. 

C 

DIMENSION FX<9)«NTX(5) 
IF( JTYPE.EQ.2) GO TO 1 
EFFsfXlLBANO) 
RETURN 
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I EFF«FXILBAND>»TEMP 
RETURN 
END 



FUNCTION WATElTENPtFXtWTXtLBANO.JTVPE) 

C 

C FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION 

C OF FREQUENCY* 

C 

OINENSION FX15)tNTXC5» 
IFCJTYPE.EG.2) 60 TO 1 
WATEsWTXCLBANDI 
RETURN 

X IF(FX(LBANDJ.LT.Q.0001» 60 TO 2 

WAT£=W TX ( LB AND l/TENP 

RETURN 
2 WATEsWTX(LBANO) 

RETURN 

ENO 



SUBROUTINE ERROR 
PRINT 1 

1 FORMAT (• •♦**•***•**• ERROR IN INPUT DATA ♦••••*•***•) 

STOP 
ENO 



SUBROUTINE REMEZCEOGEtNBANQS) 

C 

C THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM 
C FOR THE WEIGHTED CHE8YCHEV APPROXIMATION OF A CONTINUOUS 
C FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE 
C ARE A DENSE GRIO WHICH REPLACES THE FREQUENCY AXIS* THE 
C OESIREO FUNCTION ON THIS GRIO* THE WEIGHT FUNCTION ON THE 
C GRID* THE NUMBER OF COSINES* AND All INITIAL GUESS OF THE 
C EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHE8YCHEV 
C ERROR BY DETERMINING THE BEST LOCATION OF THE tXTREMAL 
C FREQUENCIES (POINTS OF MAXIMUM ERROR} ANO THEN CALCULATES 
C THE COEFFICIENTS OF THE BEST APPROXIMATION. 
C 

COMMON PI 2 1 AO • OEV « X • Y « GRIO t OES t NT t ALPHA • XEXT t NFCNS , NCR I D 
DIMENSION EDGE (20 J 

DIMENSION IEXT<661*A0<66)«ALPHA(66)*Xl£6>«T(66) 
OIMENSION DEStlQ45>«GRi0tlQ4S)tWTU043> 
DIMENSION A(66I*P165)«Q(6S) 
DOUBLE. PRECISION PI2»0NUMt DOEN.DTENP* A«P«Q 
DOUBLE PRECISION AO«DEV«X«Y 

C 

C THE PROGRAM ALLOWS A MAXIMUM NUMBER OF ITERATIONS OF 25 

C 

ITRMAXs25 
OEVL=-1«0 
N2=NFCNS*X 
NZZsNFCNS+2 
NITERsO 
100 CONTINUE 

IEXT(N2ZjtN6RID+l 
NlTERaNITER+I 

IFINITER.GT.ITRMAXI GO TO 000 

00 119 0=ltNZ 

OTENPsGRIOUEXTIJl) 

OTEMP=OCOS I 0TEMP*PI2 I 
110 X<J)-OTEMP 

JETs(NFCNS~XI/XS+X 

DO 120 UsitNZ 
120 AO< J)=Q< J«NZ«UCT> 

ONUNsO.O 

0OEN=0.0 

K=l 

DO 130 J=X*NZ 
LsXEXTU) 

OTEMP&AO(J>*DESIL) 
DNUMsONUN*OTEMP 
OTEKP=K*AOCJ)/WT(L> 
OOEN=UDEM+OTUIP 
ISO K=*K 

UEVsbNOM/OOEN 
NU-X 

IFC0EV.6T.0.Q) NU=-1 

OEV=-NU*0LV 

KzNU 

DO 1*»0 J=X«NZ 

L=X£XT|0) 

UTE«P=K*OEV/WT<L> 

YCJIcULSiD+UTE**' 
X4« Ks-K 

XFCDEV.6E.0EVL1 GO TO 190 

CALL OUCH 

GO TO 400 
150 OEVLsOEV 

JCHNGEsO 

KldEXTU) 

KNZ=IEXT<NZ> 

KLQWsQ 

MUT=-NU 

0=1 



C SEARCH FOR THE EXTREMAL FREQUENCIES OF I HE BEST 
C APPROXIMATION 

C 

200 IF(J*EQ.NZ2) YNZbCOMP 
IF(O.GE.NZZ) GO TO 300 
KUP=IEXT<J*1> 
L*IEXT4JI+1 
NUTa-NUT 

IF(J.E6.2> YlsCOMP 
COMPsOEV 

iFa.GL.KUP) GO TO 220 

ERRs£EE(L*NZ) 

ERKx(£RR.DES(L>)*WT(L> 

OTEMPxMUT*ERR*COMP 

IFCDTENP.LE.O.OI GO TO 220 

COMPsNUT»ERR 
210 L=L+1 

IF(L.GC.KUP) 60 TO 215 

ERR£6EECL*NZ> 

ERR* < ERR-DES I L > I *«T ( L > 

OTEMP=NUT*ERh-COMP 

IF(OTEMP.L£.0.0> GO TO 215 

C0MP?NUT*ERft 

GO TO 210 
215 IEXTIJ|=L-1 

J=J-»X 

KLOWsL-1 

JCHN&EsdCHN&E+X 

GO TO 200 
220 L=L-»1 
225 L=L-1 

IF(L;lE*KLOW> GO TO 256 

ERR =6 EE t L*NZ I 

ERRs t ERR-DES < L ) I *WT ( L > 

DTENPa:NUT*ERR-COMP 

IFIOTEMH.GT.Q.O) GO TO 240 

IFIJCHNGE.LE.O) 60 TO 225 

GO TO 260 
2:30 COMPsNUT*ERk 
£33 L=L-1 

iFIL.LL.KLOk) GO TO 240 

CRR=GLt»L.N2> 

EKK?(EXk-uES(L> l*WT(L> 

OTLMPsNUT*EkK-COHP 

IFCOTEMF.Lt.u.OI 66 10 240 

COrff-=NUT*tRR 

GO TO 235 
*4U KLOk=itXT(J) 

IcXTlO»=L+l . 

J=J*X 

JCHNGE=JCHMGL*1 

GO TO iCO 
250 L=UXT<J>*1 

IFtdCHN&E.GT.O) GO TO 215 
255 L=L*l 

JFU.6t.KUP) GO TO 260 

ERRsGLE(L*NZ> 

ERRs I EKR-OES I L ) > »»T I L > 

OTENPsNUT»ERH-COMP 

IF(OTEMP.LE.O.O) GO 10 255 

COMPsNUT*ERR 

GO TO 210 
260 KLOUslEXTUt 

JsJ+1 

60 TO 200 
300 IF<J.GT*NZZ> GO TO 320 

IFIK1.GT.1EXTU)) KI=1£XT<1> 

IFIKNZ.LT.IEXT(NZ>1 KNZsIEXf(NZl 

NUTlsNut 

NUTs-NU 

L=0 

KUPsKl 

COMPsYNZ* ( 1*00001) 

LUCKel 
310 L=L*1 

IF(L.GE.KUP) GO TO 315 

ERRsGEEILtNZ) 

ERRs { ERR -OES I L l> • «T ( L I 

DTEMPKNUT*ERR-COMP 

IFIOTEMP.LE.0.01 60 TO 310 

COMPsttUT«CRR 

JCNZZ 

GO TO 210 
315 LUCKS6 

GO TO 325 
320 IF<LUCK.GT«9> GO TO 350 

XF(COMP.GT.Yl) YleCONP 

Kl^IEXT (NZZ) 
325 LsNGRID+l 

KLOWsKNZ 

NUTs-NUTl 

CQMP=Yl*fl«0O001) 
330 LsL-1 

XFIL.LE.KLOW) GO TO 340 

ERRsGEE(L<NZ) 

ERR=< ERR-DES < L > J *«T < L I 

OTEMP=NUT*ERR~CQMP 
. IFIOTEMP.LE.0.01 60 TO 330 

JsNZZ 

COMP=NUT*ERR 
. LUCKsLUCK^XO 
GO TO 235 
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540 IF(LUCK«£Q«6> GO TO 370 

DO 5*9 JsltNFCNS 
3*5 IEXT<NZ2-J)=lEXTiNX-Jl 

ICXTtllcKX 

60 TO 100 
350 KK£I&XT<NZ2) 

00 360 %J=1*NFCNS 
360 IEXT(o>sIEXT<J*l> 

IEXTlNZ)=KN 

60 TO 100 
370 IF(JCHN6E.GT.0> 60 tO 100 



C CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION 

C USING THE INVERSE DISCRETE FOURIER TRAWSFORJI 

C 

tOO CONTINUE 

NNlsNFCNS-1 

FSH=1«OE-06 

GTEHP=6RI0<il 

X<NZZ»=-2»0 

CN=2*NFCN$-I 

OELFsl.O/CN 

L=l 

KKK=0 

IFIE06E<I)*EO«0*0»ANO.ED6E(2*N6ANOSI.EO*0.&> KKKsX 

IFiNfCNS,LE.3» K|$K=i 

IFCKKK.EO.l* 60 TO *0S 

UTERPsOCOSI Pl2*6RI0i 1> ) 

0MURsOC0SfPl2«6RIOl«6RIO) I 

AA=2,0/(0T£HP-0NUR) . 

88=- (OTERP+ONURl / i OTERP-ONUNf 
40S CONTINUE . 

00 430 UsitNFCNS 

FTs(J-l>*0£tF 

XT=0CGSIPI2*FT> 

IFCtCKK.EO.il 60 TO 410 

XT=(XT-BB)/AA 

FT=ARC0S(XT»/PI2 
410 XE=X(L) 

IFUT;6f,XEI 60 TO 420 

IF I < XE-XT ) *LT « FSH) 60 TO 413 

L*L*1 

60 TO 410 
41S A<«flsT«E) 

60 TO 425 
420 IFUXT-XEI.LT.FSH) 60 TO 415 

6RI0UICFT - 

AIJ)«6EE<lfNZ> 
425 CONTINUE 

IF(L,6T.1> L=L-1 
430 CONTINUE 

6RI0tl>=6TENP 

DDENSP12/CN 

DO 510 UsliftFCNS 

DTENPsO.O 

QNUfts(j-l|*D6£N 

IFCNN1*LT*11 GO TO 505 

DO 500 K=1«NR1 
500 bTEKPcDTE*P+AlK+l)»DtOSlD*U**K> 
505 DTENP62»0*OTERP*AU> 
510 ALPHA < J ) sOTCKP 

00 530 J^2*NFCNS 
bSO ALPHAlul=Z*ALPHAiJ*/LN 

ALPHA ( 1 1 =ALPhA U > /CN 

IFCKKK.tO.l? 60 TO 545 

P< 1 ) =2 . 0* ALPHA < Nl-CMS I *btf+ ALPHA ( tiHX » 

P 1 2 1 =2 • 6* A A* ALPHA I NF CNS > 

U 11 ) s ALPHA INf CNS-2 ) -ALPHA ( MFCWk ) 

00 540 J=2«K*1 

IHJ.LT.NNl) 60 TO bib 

AA=0.5*AA 

8B=0»5*88 
515 CONTINUE 

P(J+11=0*0 

00 520 KxltJ 
. AtK>=P(K| 
520 PlK>s2,0*88»MK> 

Pt2l=P(£)+AU»*2»0*AA 

JMlsJ-1 

DO 525 K=ltvWU 
523 p(K»=P(M*0<k)*AA*A<K+l> 
JP 1=4+1 

00 530 K=4*JPl 
530 P(K)=P<K|«AA*A(ft?l> 

IF< J.EW.NRl) 60 TO 540 

DO 535 K=l.u 
535 Q(K»=-A<K> 

O ( 1 1 =tt 1 1 > ♦ ALPHA t NFCN5-1- J I 
540 CONTINUE 

00 543 d=ltNFL*S 
543 ALPHA<J!=P1J) 
545 CONTINUE 

IFCNFCNS.6T.3I RETURN 

ALPHA! NFCNS+1 ) =0 .0 

ALPHA t NFCNS+2 I =0 .0 

RETURN 

END 



D0U8LE PRECISION FUNCTION U<K»N«K> 



C FUNCTION TO CALCULATE THE LA6KAN6L INTERPOLATION 

C COEFFICIENTS FOR USE IN iHt FUNCTION 6EE. 

C 

COMMON PI£.« AUtUEV*Xt V t6KID*U£S*HT « ALPHAt IEA1 tNFCNS* AIGftlO 

OIRENSION lEXTI661.A0<66>*ALPHA(66>*X(66)«Tt66I 

OlNENSlON DE3tl045)*&RID(1045}<«tTll045> 

DOUBLE PRECISION AOtOEV»X«T 

DOUBLE PRECISION 0 

DOUBLE PRECISION PI2 

0=1 »0 

O^XIKI 

00 3 LsltK 

DO 2 JsLfNiN 

IFIJ-K)l*2tl 

1 Us2»0«0*(6-X< J) i 

2 CONTINUE 

3 CONTINUE 
0=1.0/0 
RETURN 
END 



DOUBLE PRECISION FUNLUUN b££(K*NJ 

C 

C FUNCTION TO EVALUATE THt FREQUENCY RESPONSE UblNb Trtt 
C LA6KAN6E INTERPOLATION FORMULA IN THt HARTCCN f NIC FOKH 
C 

COMMON PI2 1 AOtdEV t X«Y<&KlUiUES v tiT« ALPHA *!EX1 • *FLNS,f«GR10 

0IMENS1ON IEXT(66J »AU(bbl • ALPHA (66 1 .XI fab) , Ytfco) 

DIMENSION OEStl045)f t>KlU(1045l iWTUOHb) 

DOUSLL PRECISION P.C.O.AF 

OOUbcE PRECISION P12 

DOUBLE PRLCISIO^ AU*UEViX«Y 

P=0.0 

XF^tmlUtK) 

XF=Ot0S(PI2*XF> 

0=0*0 

00 1 U=1.N 
C=XF-X(J» 
C=AO(J)/C 
U=0*C 
1 P=F*C»Y1JJ 
GEEsP/U 
RETURN 
ENU 



SUBROUTINE OUCH 
PRINT 1 - 

1 FORMAT (• ************ FAILURE lo CONVERGE ****•***»*•/ 

l*0PRObAdLE CAUSE IS MACHINE ROUNDING ERROR*/ 
2«0THE IMPULSL RESPONSE NAT bE CORRECT »/ 
3*0CHECk rflTH A FREQUENCY RESPONSE* } 

RETURN 

ENO 
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Abstract— A method of phoneme recognition of connected 
speech is described. Input to the system is assumed to consist 
of the 24 continuant phonemes in connected English speech. 
The system first categorizes each successive 20-ms segment of 
the input speech utterance as either voiced fricative, voiced 
nonfricative, unvoiced fricative or no-speech, utilizing a mea- 
sure of the relative energy balance between low and high fre- 
quencies. Next, the recognition of each 20-ms segment is per- 
formed from a distribution of axis-crossing intervals of speech 
prefHtered to emphasize each f ormant frequency range. Seg- 
mentation is performed from the results of the recognition of 
each 20-ms segment and from changes in categorization. 
Finally, the results of the recognition of each 20-ms segment 
between each pair of segmentation boundaries are combined 
and the phonemic sound occurring most frequently is printed 
out. The system has been trained for a single male speaker. 
Preliminary results for this speaker and for four 3-4-s sentences 
indicate: a correct categorization decision for about 97 per- 
cent of the input 20-ms segments, a correct recognition for 
about 78 percent of the input 20-ms segments, and an overall 
correct phoneme recognition for about 87 percent of the input 
phonemes. 



I. Introduction. 

Phoneme recognition of speech by machine has 
been a subject of increasing interest in recent years. 
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As a result, numerous techniques have been devel- 
oped and applied [1 J ^[14] . In these techniques, the 
difficulties associated with achieving phoneme recog- 
nition in total generality have forced the employment 
of constraints on the input speech utterance accept- 
able by recognition systems. Such constraints include 
a limitation on the size of the vocabulary (number of 
phonemes), a limitation on the "naturalness" of the 
utterance and a limitation on the number of speakers 
acceptable by the system. The employment of these 
three constraints, with varying degrees of restriction, 
has been universal in phoneme recognition systems. 

In the system described in this paper, the input 
speech utterance is constrained to consist of the con- 
tinuant phonemes in connected English speech. 
Hence, 24 of the possible 40 or so phonemes of En- 
glish are acceptable to the recognizer. The system 
recognizes: the eleven vowels, /i, I, e, ae, A, a, :>, u, U, 
o, 3 /; Hie four voiced fricatives, /v, 6, z, 5/; the foiir 
unvoiced fricatives, /f, 8 , s, //; the three nasals /m, n, 
17/; the two semivowels, /l, */» and the null phoneme 
(no speech). It does not presently recognize: the 
vowel glides, /e, aU, al, 0 I, iU/; the consonant glides, 
/j, co/; the affricatives, /t/, d3/; the stop consonants, 
/b, d, g, p, t, k/; or the glottal fricative /h/. The 
group of phonemes to be recognized was chosen pri- 
marily as a result of the high accuracy achieved in an 
initial study when recognizing these same phonemes 
uttered in isolation [14] . It was of interest to deter- 
mine if this high accuracy of recognition could be ac- 
complished for this same group of phonemes in con- 
tinuous speech. The resulting recognition system is 
one that vocabulary restrictions can be lessened as 
methods of recognizing the remaining phonemes are 
developed and applied. The constraint on the "natu- 
ralness" of the spoken utterance acceptable to the 
system is not made. It is assumed that no attempt is 
made to enhance recognition by other than "normal" 
enunciation or ideal noise conditions. Finally, the 
system as implemented is "trained" to accept speech 
from one talker. A suitable training procedure is 
therefore required prior to recognition. 

Four sentences containing 107 phonemes were used 
as a test of the recognition system. The system re- 
sponded correctly for about 87 percent of the pho- 
nemes. It responded incorrectly for about 4.5 per- 
cent and failed to respond for about 8.5 percent of 
the phonemes. 



